\section{Mass and Stiffness Matrix in FEM}

Here is just a very simple and intuitive derivative for the formulation of mass and stiffness matrix. Actually there is no special to these two matrices in FEM. They just come out naturally when forming the linear system of \textit{Galerkin Projection} method, like many other coefficient matrices. Their names come from classical mechanics and elastic theory. We are going to show it in two very simple cases in 1D but the idea can be generalized to higher dimesion.

First, we gonna to show a function approximation case for mass matrix. Given the following equation,

\begin{equation}\nonumber
u = f
\end{equation}

where u is the solution we wanna seek to satisfy the above equation while f is some known function. Apparently the solution is f itself. However, we may not get it that easy if the left hand side becomes way more complicated, like arbitary nonliner differential algebraic equation. Thus we try to use some linear combination of a set of simple basis equation to approximate the true solution, like

\begin{equation}\nonumber
u \approx u_h = \sum_{i}c_i\phi_i
\end{equation}

The error between the true solution and the approximation is 

\begin{equation}\nonumber
e = u - u_h
\end{equation}

In order to determine the unknown coefficients $c_i$, we can convert this problem to a variational optimization problem like minimize certain norm of the error. Such procedure in this case coincide with the idea of \textit{Galerkin Projection} method which states that the error should be orthogonal with the function space \textbf{V} spanned by some test functions $\omega_i$. Here we just assume the test functions and the basis functions or trial functions are the same. Then we have the following

\begin{equation}\nonumber
<e, v> = 0, \forall v \in \textbf{V}
\end{equation}

Notice $<\cdot,\cdot>$ is a general bilinear form here related to the definition of orthogonallity, like inner product for example. Then we have

\begin{equation}\nonumber
\begin{split}
&<e, \phi_i> = 0\\
&\Rightarrow <f-\sum_{j}c_j\phi_j, \phi_i>=0 \Rightarrow <\sum_{j}c_j\phi_j,\phi_i> = <f, \phi_i>\\
&\Rightarrow \sum_{j}<\phi_j,\phi_i>c_j = <f, \phi_i>
\end{split}
\end{equation}

This is simply a linear system and $<\phi_j,\phi_i>$ forms the mass matrix. If we use $l_2$ norm or inner product (for real field only) to define the bilinear form, we have

\begin{equation}\nonumber
<\phi_j,\phi_i> = \int_{\Omega}\phi_j\phi_idx
\end{equation}

Normally the mass matrix is dense, but if choosing the basis functions and quadrature rule carefully like the combination of P1 and  Trapezoidal integration scheme, we can diagonalize or lump the mass matrix, meaning

\begin{equation}\nonumber
<\phi_j,\phi_i> = \delta_{ij}
\end{equation}

For stiffness matrix, let's take a look at one simple second order ordinary differential equation

\begin{equation}\nonumber
\frac{d^2u}{dx^2} = f
\end{equation}

Similar to the first example, we try to find a linear combination of some basis functions to approximate the true solution. However, in this case, we don't know what the true solution is, the only thing we can messure is the so called residual which is defined as follow

\begin{equation}\nonumber
r = f - \frac{d^2u_h}{dx^2}
\end{equation}

Follow the same procedure of \textit{Galerkin Project} method, we have

\begin{equation}\nonumber
\begin{split}
&<r, \phi_i> = 0\\
&\Rightarrow <f-\frac{d^2u_h}{dx^2}, \phi_i>=0\\
&\Rightarrow \sum_{j}<\frac{d^2\phi_j}{dx^2},\phi_i>c_j = <f, \phi_i>
\end{split}
\end{equation}

Basically $<\frac{d^2\phi_j}{dx^2},\phi_i>$ forms the stiffness matrix, but usually people like to write it in a different way by applying rule of integration by part. Here we use the $l_2$ norm as the bilinear operator

\begin{equation}\nonumber
\begin{split}
<\frac{d^2\phi_j}{dx^2},\phi_i> &= \int_{\Omega}\frac{d^2\phi_j}{dx^2}\phi_idx = -\int_{\Omega}\frac{d\phi_j}{dx}\frac{d\phi_i}{dx}dx + \left.[\frac{d\phi_j}{dx}\phi_i]\right\vert_{\Gamma}\\
&= -\int_{\Omega}\frac{d\phi_j}{dx}\frac{d\phi_i}{dx}dx = - <\frac{d\phi_j}{dx},\frac{d\phi_i}{dx}>
\end{split}
\end{equation}